# Turn on the gmri font for plots
showtext::showtext_auto()
# File paths for various extents based on params$region
region_paths <- get_timeseries_paths(region_group = "gmri_sst_focal_areas", 
                                     mac_os = "mojave")

# Load the bounding box for Andy's GOM to show they align
poly_path     <- region_paths[[params$region]][["shape_path"]]
region_extent <- st_read(poly_path, quiet = TRUE)


# Load the timeseries
timeseries_path <- region_paths[[params$region]][["timeseries_path"]]
region_timeseries <- read_csv(timeseries_path, 
                              col_types = cols(), 
                              guess_max = 1e6) %>% 
  mutate(tod = format(time, format = "%H:%M:%S")) # %>% 
  # filter(tod != "12:00:00") # removes the NCEI dates that leak through
# Clean up the data - add labels
region_timeseries <- region_timeseries %>% 
  mutate(time = as.Date(time)) %>% 
  distinct(time, .keep_all = T) %>% 
  supplement_season_info() %>% 
  filter(year %in% c(1982:2022))



# Summarize by year to return mean annual anomalies and variance
annual_summary <- region_timeseries %>% 
  filter(year %in% c(1982:2021)) %>% 
  group_by(year) %>% 
  temperature_summaries() %>% 
  mutate(yr_as_dtime = as.Date(paste0(year, "-07-02")),
         anom_direction = ifelse(area_wtd_anom > 0, "Hot", "Cold"))


# Get heatwave statuses for each day:
# Uses area weighted sst by default
region_hw <- pull_heatwave_events(
  temperature_timeseries = region_timeseries, 
  threshold = 90, 
  clim_ref_period = c("1982-01-01", "2011-12-31")) %>% 
  supplement_hw_data() %>% 
  filter(doy != 366)



 # Last Calendar Year
this_yr <- region_hw %>% 
  filter(year(time) == 2022)


# Set up the date within a year for Heatmap
base_date <- as.Date("2000-01-01")
grid_data <- region_hw %>% 
  mutate(year = year(time),
         yday = yday(time),
         flat_date = as.Date(yday-1, origin = base_date),
         yr_rev = factor(year),
         yr_rev = fct_rev(yr_rev))


# # Global Temperature Anomaly Rates
global_anoms <- read_csv(
    paste0(oisst_path, "global_timeseries/global_anoms_1982to2011.csv"), 
    guess_max = 1e6,
    col_types = cols()) %>% 
  mutate(year = year(time)) %>% 
  filter(between(year, 1982, 2021))

# summarize by year again
global_summary <- global_anoms %>% 
  group_by(year) %>% 
  temperature_summaries() %>% 
  mutate(yr_as_dtime = as.Date(paste0(year, "-07-02")),
         anom_direction = ifelse(area_wtd_anom > 0, "Hot", "Cold"))


####  Run Warming Rate regressions  ####
rate_data <- list(
  "GoM All F" = get_decadal_rates(temp_df = annual_summary, 
                                temp_col = "area_wtd_f", 
                                year_col = "year", 
                                year_lim = c(1982, 2021), 
                                area_name = "GoM", 
                                degree_c = F),
  "GoM All C" = get_decadal_rates(temp_df = annual_summary, 
                                temp_col = "area_wtd_sst", 
                                year_col = "year", 
                                year_lim = c(1982, 2021), 
                                area_name = "GoM", 
                                degree_c = T),
  "GoM 15" = get_decadal_rates(temp_df = annual_summary, 
                                temp_col = "area_wtd_f", 
                                year_col = "year", 
                                year_lim = c(2007, 2021), 
                                area_name = "GoM", 
                                degree_c = F),
  "Global All" = get_decadal_rates(temp_df = global_summary, 
                                temp_col = "area_wtd_f", 
                                year_col = "year", 
                                year_lim = c(1982, 2021), 
                                area_name = "Global", 
                                degree_c = F),
  "Global All C" = get_decadal_rates(temp_df = global_summary, 
                                temp_col = "area_wtd_sst", 
                                year_col = "year", 
                                year_lim = c(1982, 2021), 
                                area_name = "Global", 
                                degree_c = T))

Gulf Of Maine Sea Surface Temperature Outlook

This report was created to track the sea surface temperature regimes for marine regions of interest to the Gulf of Maine Research Institute. The default region being a central snapshot of the Gulf of Maine.

Satellite sea surface temperature data used was obtained from the National Center for Environmental Information (NCEI). With all maps and figures displaying NOAA’s Optimum Interpolation Sea Surface Temperature Data.

Region Extent

The spatial extent for Gulf Of Maine is displayed below. This bounding box is the same bounding box coordinates used to clip the OISST data when constructing the time series data from the array.

# Pull extents for the region to set crop extent
crop_x <- st_bbox(region_extent)[c(1,3)]
crop_y <- st_bbox(region_extent)[c(2,4)]


# Map the study area
map_study_area(region_extent = region_extent,
               x_buffer =  c(-2, 2),
               y_buffer =  c(-0.75, 0.75), 
               new_england_sf = new_england, 
               canada_sf = canada,
               greenland_sf = greenland)

Regional Timeseries

Area-specific time series are the most basic building block for relaying temporal trends. For any desired area (represented by a spatial polygon) we can generate a time series table of the mean sea surface temperature within that area for each day. Additionally, we can compare how observed temperatures correspond with the expected conditions based on a climatology using a specified reference period.

# Display Table of last 6 entries
tail(region_timeseries) %>% 
  mutate(across(where(is.numeric), round, 2)) %>% 
  select(
    Date = time,
    `Sea Surface Temperature` = sst,
    #`Area-Weighted SST` = area_wtd_sst,
    `Day of Year` = modified_ordinal_day,
    `Climate Avg.` = sst_clim,
    #`Area-weighted Climate` = area_wtd_clim,
    `Temperature Anomaly` = sst_anom#,
    #`Area-Weighted Anomaly` = area_wtd_anom
    
  ) %>% gt() %>% 
    tab_header(
    title = md(paste0("**", tidy_name, " - Regional Sea Surface Temperature", "**")), 
    subtitle = paste("Temperature Unit: Celsius")) %>%
  tab_source_note(
    source_note = md("*Data Source: NOAA OISSTv2 Daily Sea Surface Temperature Data.*") ) %>% 
  tab_source_note(md("*Climatology Reference Period: 1982-2011.*"))
Gulf Of Maine - Regional Sea Surface Temperature
Temperature Unit: Celsius
Date Sea Surface Temperature Day of Year Climate Avg. Temperature Anomaly
2022-05-19 10.32 140 7.86 2.47
2022-05-20 10.03 141 7.98 2.06
2022-05-21 10.36 142 8.17 2.19
2022-05-22 10.57 143 8.32 2.25
2022-05-23 10.80 144 8.46 2.34
2022-05-24 10.55 145 8.63 1.92
Data Source: NOAA OISSTv2 Daily Sea Surface Temperature Data.
Climatology Reference Period: 1982-2011. 
# march 1st sst
mar1 <- region_timeseries %>% 
  filter(modified_ordinal_day == 61) %>% 
  distinct(sst_clim) %>% 
  pull(sst_clim)

Each of our Climatologies are currently set up to calculate daily averages on a modified year day, such that every March 1st and all days after fall on the same day, regardless of whether it is a leap year or not.

This preserves comparisons across calendar dates such-as: “The average temperature on march 1st is 4.2378197` for the reference period 1982 to 2011”

In these tables Sea Surface Temperature is the mean temperature observed for that date averaged across all cells within the area. Climate Avg. & Climate SD are the climate means and standard deviations for a 1982-2011 climatology. Temperature Anomaly is the daily observed sea surface temperature minus the climate mean.

Warming Rates

Regional warming trends below were calculated using all the available data for complete years beginning with 1982 through the end of 2021. The overlaid trend lines then track how warming has increased with time. A dotted line has been included to show how the global average temperature has changed during the same period.

Annual

# Build Regression Equation Labels
eq_all_c  <- rate_data[["GoM All C"]][["eq_label"]]
eq_all    <- rate_data[["GoM All F"]][["eq_label"]]
eq_15     <- rate_data[["GoM 15"]][["eq_label"]]
eq_global <- rate_data[["Global All"]][["eq_label"]]
eq_global_c <- rate_data[["Global All C"]][["eq_label"]]



# Generate a smoothed temperature line using splines
yearly_temp_smooth <-  as.data.frame(spline(annual_summary$yr_as_dtime, annual_summary$area_wtd_anom)) %>% 
  mutate(x = as.Date(x, origin = "1970-01-01"))
####  Annual Trend Plot  ####
ggplot(data = annual_summary, aes(yr_as_dtime, area_wtd_anom)) +
  
  # Add daily data
  geom_line(data = region_timeseries,
            aes(time, area_wtd_anom, color = "Daily Temperatures")) +
  
  # Overlay yearly means
  geom_line(data = yearly_temp_smooth, aes(x, y, color = "Average Yearly Temperature"), alpha = 0.7, linetype = 2) +
  # geom_line(color = "gray10", size = 1) +
  geom_point(color = "gray10", alpha = 0.7, size = 0.75) +
  
  # Add regression lines
  geom_textsmooth(data = filter(global_summary, year <= 2021),
              method = "lm", text_smoothing = 30,
              label = "Global Trend",
              color = gmri_cols("green"),
              linewidth = 1,
              formula = y ~ x, se = F,
              linetype = 3, hjust = 0.925) +
 
  geom_smooth(data = filter(annual_summary, year <= 2021),
              method = "lm",
              aes(color = "1982-2021 Regional Trend"), #label = "40-Year Trend",
              formula = y ~ x, se = F,
              linetype = 1) +
  # geom_smooth(data = filter(annual_summary, year %in% c(2007:2021)),
  #             
  #             method = "lm", 
  #             aes(color = "2007-2021 Regional Trend"),
  #             formula = y ~ x, se = F, 
  #             linetype = 2) +

  
  # Manually add equations so they show yearly not daily coeff
  geom_text(data = data.frame(), 
            aes(label = eq_all, x = min(region_timeseries$time), y = Inf),
            hjust = 0, vjust = 3, color = gmri_cols("gmri blue")) +
  geom_text(data = data.frame(), 
            aes(label = eq_global, x = min(region_timeseries$time), y = Inf),
            hjust = 0, vjust = 4.5, color = gmri_cols("green")) +
  # geom_text(data = data.frame(), 
  #           aes(label = eq_15, x = min(region_timeseries$time), y = Inf),
  #           hjust = 0, vjust = 6, color = gmri_cols("orange")) +

  # Colors
  scale_color_manual(values = c(
    "1982-2021 Regional Trend" = as.character(gmri_cols("gmri blue")),
    #"2007-2021 Regional Trend" = as.character(gmri_cols("orange")),
    #"1982-2021 Global Trend"   = as.character(gmri_cols("green")),
    "Average Yearly Temperature" = "gray10",
    "Daily Temperatures"       = "gray90")) +
    
  # Axes
   scale_y_continuous(
     sec.axis = sec_axis(trans = ~as_fahrenheit(., data_type = "anomalies"),
                         labels =  number_format(suffix = " \u00b0F")),
     labels = number_format(suffix = " \u00b0C")) +
  # theme
   theme(legend.title = element_blank(),
        legend.position = "bottom",
        legend.background = element_rect(fill = "transparent"),
        legend.key = element_rect(fill = "transparent", color = "transparent"),
        #panel.grid = element_blank()
        ) +
  # labels + theme
  labs(x = "", 
       y = "Sea Surface Temperature Anomaly",
       caption = paste0("Anomalies calculated using 1982-2011 reference period."))

Without Daily Temperatures

####  Annual Trend Plot  ####
ggplot(data = annual_summary, aes(yr_as_dtime, area_wtd_anom)) +
  
  # Overlay yearly means and lines connecting them
  geom_col(aes(fill = "Yearly Anomaly"), alpha = 0.7, size = 0.75) +
  scale_fill_manual(values = c("Yearly Anomaly" = "gray90")) +
  
  # Add regression lines
  geom_textsmooth(data = filter(global_summary, year <= 2021),
              method = "lm", text_smoothing = 30,
              size = 4,
              label = "Global Trend",
              color = gmri_cols("green"),
              linewidth = 1,
              formula = y ~ x, se = F,
              linetype = 1, hjust = 0.925) +
 
  geom_textsmooth(data = filter(annual_summary, year <= 2021),
              method = "lm",
              # aes(color = "1982-2021 Regional Trend"), #label = "40-Year Trend",
              size = 4, linewidth = 1,
              label = "1982-2021 Regional Trend",
              formula = y ~ x, se = F,
              linetype = 1, linewidth = 1, color = gmri_cols("gmri blue"),
              hjust = 0.50, vjust = -1.2
              ) +
  geom_textsmooth(data = filter(annual_summary, year %in% c(2007:2021)),
              method = "lm",
              # aes(color = "1982-2021 Regional Trend"), #label = "40-Year Trend",
              size = 4, linewidth = 1,
              label = "2007-2021 Regional Trend",
              formula = y ~ x, se = F,
              linetype = 1, linewidth = 1, color = gmri_cols("orange"),
              #boxlinetype = "dotted", boxlinewidth = 0.5,
              hjust = 0.95, vjust = -1.2
              ) +
   
  
  # Manually add equations so they show yearly not daily coeff
  geom_text(data = data.frame(), 
            aes(label = eq_all, x = min(region_timeseries$time), y = Inf),
            hjust = 0, vjust = 6, color = gmri_cols("gmri blue")) +
  geom_text(data = data.frame(), 
            aes(label = eq_15, x = min(region_timeseries$time), y = Inf),
            hjust = 0, vjust = 7.5, color = gmri_cols("orange")) +
  geom_text(data = data.frame(), 
            aes(label = eq_global, x = min(region_timeseries$time), y = Inf),
            hjust = 0, vjust = 9, color = gmri_cols("green")) +
    
  # Axes
   scale_y_continuous(sec.axis = sec_axis(trans = ~as_fahrenheit(., data_type = "anomalies"),
                                         labels =  number_format(suffix = " \u00b0F")),
                     labels = number_format(suffix = " \u00b0C")) +
  
  
  # labels + theme
  labs(title = "Gulf of Maine Warming Faster then Global Average",
       x = "", 
       y = "Sea Surface Temperature Anomaly",
       caption = paste0("Anomalies calculated using 1982-2011 reference period.")) +
  # theme
   theme(legend.title = element_blank(),
        legend.position = c(0.85, 0.1),
        legend.background = element_rect(fill = "transparent"),
        legend.key = element_rect(fill = "transparent", color = "transparent"))

Quarterly

# Doing seasons by meteorological Definitions
quarter_summary <- region_timeseries %>% 
 filter(year >= 1982) %>% 
  group_by(year = season_yr, season_eng) %>% 
  summarise(sst = mean(sst, na.rm = T),
            sst_anom = mean(sst_anom, na.rm = T), 
            area_wtd_sst = mean(area_wtd_sst, na.rm = T),
            area_wtd_anom = mean(area_wtd_anom, na.rm = T),
            .groups = "drop") %>% 
  mutate(season_eng = factor(season_eng, levels = c("Winter", "Spring", "Summer", "Fall")))

# Plot
quarter_summary %>% 
  ggplot(aes(year, area_wtd_anom)) +
  geom_line(group = 1, color = "gray60", linetype = 1) +
  geom_point(size = 0.75) +
  geom_smooth(method = "lm", 
              aes(color = "Regional Trend"),
              formula = y ~ x, se = F, linetype = 2) +
  stat_poly_eq(formula = y ~ x,
               color = gmri_cols("orange"),
               aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
               parse = T) +
  scale_color_manual(values = c("Regional Trend" = as.character(gmri_cols("orange")))) +
  labs(x = "", 
       y = "Sea Surface Temperature Anomaly",
       caption = "Regression coefficients reflect annual change in sea surface temperature.") +
   scale_y_continuous(sec.axis = sec_axis(trans = ~as_fahrenheit(., data_type = "anomalies"),
                                         labels =  number_format(suffix = " \u00b0F")),
                     labels = number_format(suffix = " \u00b0C")) +
  theme(legend.title = element_blank(),
        legend.position = "none",
        legend.background = element_rect(fill = "transparent"),
        legend.key = element_rect(fill = "transparent", color = "transparent")) +
  facet_wrap(~season_eng, ncol = 1)

Overall Temperature Increase

dat_region <- annual_summary %>% 
  filter(year %in% c(1982, 2021)) 
dat_global <- global_summary %>% 
  filter(year %in% c(1982, 2021)) 
dat_list <- list(dat_region, dat_global) %>% setNames(c(tidy_name, "Global Oceans"))
dat_combined <- bind_rows(dat_list, .id = "Area") %>% 
  select(Area, area_wtd_sst, year) %>% 
  pivot_wider(names_from = year, values_from = area_wtd_sst)

ggplot(dat_combined, aes(x = `1982`, xend = `2021`, y = fct_rev(Area))) +
  geom_dumbbell(colour = "lightblue", 
                colour_xend = gmri_cols("gmri blue"), 
                size = 3, 
                dot_guide = TRUE, 
                dot_guide_size = 0.5) +
  geom_text_repel(
    aes(x = `2021`, y = fct_rev(Area), 
        label = str_c("+", round(`2021`- `1982`, 2), " C")),
                color = gmri_cols("gmri blue"),
                vjust = 4,
                hjust = 0,
                segment.size = 0.5,
                seed = 123) +
  labs(x = "Sea Surface Temperature", 
       title = "Change in Sea Surface Temeprature - 1982-2021", 
       y = "") +
   scale_x_continuous(sec.axis = sec_axis(trans = ~as_fahrenheit(., data_type = "temperature"),
                                         labels =  number_format(suffix = " \u00b0F")),
                     labels = number_format(suffix = " \u00b0C"))

Long-Term Trend

Thanks to some long-running temperature monitoring programs we are able to track monthly sea surface temperatures as far back as 1850. The figure below shows how the region’s temperatures have fluctuated as measured by the Extended Reconstructed Sea Surface Temperature (ERSST) record.

Once we hit 1982 the availability of a higher resolution temperature resource becomes available in the form of the Optimum Interpolation Sea Surface Temperature (OISST) record. The ERSST data is of a coarser resolution (monthly measurements at a 0.5 x 0.5 degree grid) and does not capture the warmer inshore dynamics (a bias to show colder temperature), but is within half of a degree of the more modern measurements.

# # Pathto ERSST on box
# Updated 2022-02-4
ersst_path <- cs_path("res", "ERSSTv5")
ersst_tl <- read_csv(str_c(ersst_path, "ERSSTv5_anom_apershing_gulf_of_maine.csv")) 

# format
ersst_tl <- ersst_tl %>% 
  rename(Date = time) %>% 
  mutate(yr = year(Date),
         sst_clim = sst - sst_anom,
         area_wtd_clim = area_wtd_sst - area_wtd_sst_anom)


# Get Yearly Average
ersst_yr <- ersst_tl %>% 
  group_by(yr) %>% 
  summarise(sst = mean(sst, na.rm = T), 
            sst_anom = mean(sst_anom, na.rm = T),
            .groups = "drop") %>% 
  mutate(yr = as.numeric(yr)) %>% 
  filter(yr <= 2021) 


# Same for OISST
oisst_yr <- region_hw %>% 
  filter(between(yr, 1982, 2021)) %>% 
  group_by(yr) %>% 
  summarise(sst = mean(sst), .groups = "drop")

# Make Splines
ersst_smooth <- as.data.frame(spline(ersst_yr$yr, ersst_yr$sst))
oisst_smooth <- as.data.frame(spline(oisst_yr$yr, oisst_yr$sst))


# Plot them
ggplot() +
  geom_line(data = ersst_smooth, aes(x, y, color = "ERSSTv5"),
            group = 1, alpha = 0.4, linetype = 1) +
  geom_point(data = ersst_yr, aes(yr, sst, color = "ERSSTv5"), 
             size = 1) +
 geom_line(data = oisst_smooth, aes(x, y, color = "OISSTv2"),
            group = 1, alpha = 0.4, linetype = 1) +
  geom_point(data = oisst_yr, aes(yr, sst, color = "OISSTv2"), 
             size = 1) +
  scale_color_gmri() +
  theme_gmri() +
  labs(x = "", y = "Sea Surface Temperature", 
       color = "Temperature Data Record",
       title = "Long-Term Temperature Record for the Gulf of Maine",
       subtitle = "Current temperatures above 150-year highs") +
  scale_y_continuous(sec.axis = sec_axis(trans = ~as_fahrenheit(., data_type = "anomalies"),
                                         labels =  number_format(suffix = " \u00b0F")),
                     labels = number_format(suffix = " \u00b0C")) +
  guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) +
  facet_zoom(xlim = c(1982, 2021)) +
  theme(legend.position = "bottom")

Long term temperature record for the Gulf of Maine

Marine Heatwaves

For the figures below heatwave events were determined using the methods of Hobday et al. 2016 and implemented using the R package {heatwaveR}.

A marine heatwave is defined as a situation when seawater temperatures exceeds a seasonally-varying threshold (usually the 90th percentile) for at least 5 consecutive days. Successive heatwaves with gaps of 2 days or less are considered part of the same event. The heatwave threshold used below was 90%. The heatwave history for Gulf Of Maine is displayed below:

Interactive SST Timeline

For anything we wish to host on the web there is an option to display tables and graphs that are interactive. Interactivity allows users to pan, zoom, and highlight discrete observations.

# # Option 1: Grab last 365 days
# 
# # Grab data from the most recent year through present day to plot
# last_year <- Sys.Date() - 365 #1 year from current date
# last_yr_heatwaves <- region_hw %>% filter(time >= last_year)
# last_yr <- year(last_yr)
# # Option 2: Grab last Full Year:
# last_year <- Sys.Date() - 365 #1 year from current date
# 
# # wind it back to first day of the year
# last_year <- last_year - yday(last_year) + 1 
# 
# # Filter out:
# last_yr_heatwaves <- region_hw %>% filter(year(time) == year(last_year))
# last_yr <- year(last_yr)
# Option 3: Last complete year
last_year <- 2021
last_yr_heatwaves <-  region_hw %>% filter(year(time) == last_year)
# Get number of heatwave events and total heatwave days for last year
# How many heatwave events:

# How many heatwave events
num_events   <- max(last_yr_heatwaves$mhw_event_no, na.rm = T) - min(last_yr_heatwaves$mhw_event_no, na.rm = T) + 1

# Data for only the current year
this_yr_hw <- region_hw %>% filter(year(time) == year(Sys.Date()))

# Number of heatwave events
num_hw_days <- sum(this_yr_hw$mhw_event, na.rm = T)
# Plot the interactive timeseries
last_yr_heatwaves %>% 
  plotly_mhw_plots()

2022 Temperatures

# Make Plot
hw_temp_p <- year_hw_temps(year_hw_dat = this_yr_hw, temp_units = "C")

# Show Figure
hw_temp_p

2022 Anomalies

# Show figure
hw_anom_p <- year_hw_anoms(year_hw_dat = this_yr_hw, temp_units = "C")
hw_anom_p

Heatwave Events

# Plot heatmap
heatwave_heatmap <- heatwave_heatmap_plot(hw_dat = grid_data, temp_units = "C")


# Assemble pieces
heatwave_heatmap

SST Anomaly Maps

The 2021 global sea surface temperature anomalies have been loaded and displayed below to visualize how different areas of the ocean experience swings in temperature.

# Access information to netcdf on box
nc_year <- "2021"
anom_path <- str_c(oisst_path, "annual_anomalies/1982to2011_climatology/daily_anoms_", nc_year, ".nc")

# Load 2021 as stack
anoms_2021 <- stack(anom_path)

Global Anomalies

Annual Map

# Make Annual Average
ann_anom_ras <- calc(anoms_2021, fun = mean, na.rm = T)

# Color scale title
sst_lab <- "Sea Surface Temperature Anomaly"

# Color limit for palettes
temp_limits <- c(-2, 2)
temp_breaks <- c(temp_limits[1], temp_limits[1]/2,  0, temp_limits[2]/2, temp_limits[2])
temp_labels <- str_c(c(str_c("< ", temp_limits[1]), temp_limits[1]/2, 0, temp_limits[2]/2, str_c("> ", temp_limits[2])), "\u00b0C")

# Make st object
ann_anom_st <- st_as_stars(rotate(ann_anom_ras)) 

# Build Map
ggplot() +
  geom_stars(data = ann_anom_st) +
  geom_sf(data = world, fill = "gray30", color = "white", size = .25) +
  scale_fill_distiller(
    palette = "RdBu", 
    na.value = "transparent", 
    limit = temp_limits, 
    oob = scales::squish,
    breaks = temp_breaks, 
    labels = temp_labels) +
  guides("fill" = guide_colorbar(
    title = sst_lab, 
    title.position = "right", 
    title.hjust = 0.5,
    barheight = unit(3, "inches"), 
    frame.colour = "black", 
    ticks.colour = "black")) +  
  coord_sf(expand = F) +
  map_theme +
  theme(legend.position = "right", 
        legend.title = element_text(angle = 90)) +
  labs(title = str_c("Global Temperature Anomalies - ", nc_year))

Seasonal Maps

#### Making Averages for Seasons

# Get previous year winter
last_nc_yr <- as.numeric(nc_year) - 1
last_yr_anom_path <- str_c(oisst_path, "annual_anomalies/1982to2011_climatology/daily_anoms_", last_nc_yr, ".nc")

# Load 2020 as stack to get the full winter
last_yr_anoms <- stack(last_yr_anom_path)

# Join to 2021
anoms_double <- stack(list(last_yr_anoms, anoms_2021))


# Drop december 2021 so its not influencing the Winter of 2020-2021
dec_2021 <- "X2021.12"
not_dec21 <- which(str_sub(names(anoms_double), 1, 8) != dec_2021)
anoms_nodec <- anoms_double[[not_dec21]]


# Set up list to use map()
month_key <- list("Winter" = c("12", "01", "02"),
                  "Spring" = c("03", "04", "05"),
                  "Summer" = c("06", "07", "08"),
                  "Fall"   = c("09", "10", "11"))



#  Get mean anoms across seasons
season_stacks <- map(month_key, function(mon){
  
  # Get names from the stack
  stack_names <- names(anoms_nodec)
  stack_months <- str_sub(stack_names, 7,8)
  
  # layers with correct months:
  in_season <- which(stack_months %in% mon)
  
  # Get mean across those months
  season_mean <- calc(anoms_nodec[[in_season]], mean, na.rm = T)
  # season_mean <- st_as_stars(rotate(season_mean))
  return(season_mean)
  
  
})



# Plotting the seasons
seas_anom_maps <- function(season, lab_years, temp_lim = 5){
  
  #temp_limits <- max(abs(values(season_stacks[[season]])), na.rm = T) * c(-1,1) 
  
  # Center the color scale with Dynamic Limits
  temp_limits <- c(-1,1) * temp_lim
  temp_breaks <- c(temp_limits[1], temp_limits[1]/2,  0, temp_limits[2]/2, temp_limits[2])
  temp_labels <- str_c(c(str_c("< ", temp_limits[1]), temp_limits[1]/2, 0, temp_limits[2]/2, str_c("> ", temp_limits[2])), "\u00b0C")
  
  # Make map
  seas_map <- ggplot() +
    geom_stars(data = st_as_stars(rotate(season_stacks[[season]]))) +
    geom_sf(data = world, fill = "gray30", color = "white", size = .25) +
    scale_fill_distiller(palette = "RdBu", 
                     na.value = "transparent", 
                     limit = temp_limits, 
                     oob = scales::squish,
                     breaks = temp_breaks, 
                     labels = temp_labels) +
    guides("fill" = guide_colorbar(title = sst_lab, 
                                   title.position = "right", 
                                   title.hjust = 0.5,
                                   barheight = unit(3, "inches"), 
                                   frame.colour = "black", 
                                   ticks.colour = "black")) +  
    coord_sf(expand = F) +
   map_theme +
    theme(legend.position = "right", 
          legend.title = element_text(angle = 90)) +
    labs(title = str_c(season, " - ", lab_years))
  return(seas_map)

  
}
Winter
seas_anom_maps("Winter", "2020-2021", temp_lim = 2)

Spring
seas_anom_maps("Spring", "2021", temp_lim = 2)

Summer
seas_anom_maps("Summer", "2021", temp_lim = 2)

Fall
seas_anom_maps("Fall", "2021", temp_lim = 2)

Regional Maps

# # Expand the area out to see the larger patterns
# crop_x <- crop_x + c(-2.5, 2.5)
# crop_y <- crop_y + c(-1.5, 1)

# Make a new extent for cropping
region_extent_expanded <- st_sfc(st_polygon(list(
  rbind(c(crop_x[[1]], crop_y[[1]]),  
        c(crop_x[[1]], crop_y[[2]]), 
        c(crop_x[[2]], crop_y[[2]]), 
        c(crop_x[[2]], crop_y[[1]]), 
        c(crop_x[[1]], crop_y[[1]])))))

region_extent_expanded <- st_as_sf(region_extent_expanded)

Annual Map

# Mask the annual average
reg_ann_anom <- mask_nc(ras_obj = ann_anom_ras, mask_shape = region_extent_expanded)

# Center temperature limits
# temp_limits <- max(abs(values(reg_ann_anom)), na.rm = T) * c(-1,1) # Dynamic Limits

# Color limit for palettes
temp_limits <- c(-5, 5)
temp_breaks <- c(temp_limits[1], temp_limits[1]/2,  0, temp_limits[2]/2, temp_limits[2])
temp_labels <- str_c(c(str_c("< ", temp_limits[1]), temp_limits[1]/2, 0, temp_limits[2]/2, str_c("> ", temp_limits[2])), "\u00b0C")

# Build Map
ggplot() +
  geom_stars(data = st_as_stars(reg_ann_anom)) +
  geom_sf(data = new_england, fill = "gray30", color = "white", size = .25) +
  geom_sf(data = canada, fill = "gray30", color = "white", size = .25) +
  geom_sf(data = greenland, fill = "gray30", color = "white", size = .25) +
  scale_fill_distiller(palette = "RdBu", 
                     na.value = "transparent", 
                     limit = temp_limits, 
                     oob = scales::squish,
                     breaks = temp_breaks, 
                     labels = temp_labels) +
    guides("fill" = guide_colorbar(title = sst_lab, 
                                   title.position = "right", 
                                   title.hjust = 0.5,
                                   barheight = unit(3, "inches"), 
                                   frame.colour = "black", 
                                   ticks.colour = "black")) +  
  map_theme +
  theme(legend.position = "right", 
        legend.title = element_text(angle = 90)) +
  coord_sf(xlim = crop_x, 
           ylim = crop_y, expand = F) +
  labs(title = str_c("Regional Temperature Anomalies - ", nc_year))

Seasonal Maps

# Mask the Seasons
seasons_masked <- map(season_stacks, mask_nc, region_extent_expanded)

# Plot the seasons
plot_masked_season <- function(season, year_lab, temp_lim = 5){
  
  # Grab Season
  reg_seas_anom <- seasons_masked[[season]]
  
  # Get temp limits
  # temp_limits <- max(abs(values(reg_seas_anom)), na.rm = T) * c(-1,1) # Dynamic Limits
  # Color limit for palettes
  temp_limits <- c(-1, 1) * temp_lim
  temp_breaks <- c(temp_limits[1], temp_limits[1]/2,  0, temp_limits[2]/2, temp_limits[2])
  temp_labels <- str_c(c(str_c("< ", temp_limits[1]), temp_limits[1]/2, 0, temp_limits[2]/2, str_c("> ", temp_limits[2])), "\u00b0C")
  
  
  # Build Map
  seas_map <- ggplot() +
    geom_stars(data = st_as_stars(reg_seas_anom)) +
    geom_sf(data = new_england, fill = "gray30", color = "white", size = .25) +
    geom_sf(data = canada, fill = "gray30", color = "white", size = .25) +
    geom_sf(data = greenland, fill = "gray30", color = "white", size = .25) +
    scale_fill_distiller(palette = "RdBu", 
                     na.value = "transparent", 
                     limit = temp_limits, 
                     oob = scales::squish,
                     breaks = temp_breaks, 
                     labels = temp_labels) +
    guides("fill" = guide_colorbar(title = sst_lab, 
                                   title.position = "right", 
                                   title.hjust = 0.5,
                                   barheight = unit(3, "inches"), 
                                   frame.colour = "black", 
                                   ticks.colour = "black")) +  
  map_theme +
  theme(legend.position = "right", 
        legend.title = element_text(angle = 90)) +
    coord_sf(xlim = crop_x, 
             ylim = crop_y, expand = F) +
    labs(title = str_c(season, " - ", year_lab))
  
  return(seas_map)
  
}
Winter
plot_masked_season("Winter", year_lab = "2020-2021", temp_lim = 3.5)

Spring
plot_masked_season("Spring", year_lab = "2021", temp_lim = 3.5)

Summer
plot_masked_season("Summer", year_lab = "2021", temp_lim = 3.5)

Fall
plot_masked_season("Fall", year_lab = "2021", temp_lim = 3.5)

Heatwave Progression

Looking specifically at the last heatwave event, we can step through how the event progressed over time, and developing pockets or warmer/colder water masses.

# Identify the last heatwave event that happened
last_event <- max(region_hw$mhw_event_no, na.rm = T)

# Pull the dates of the most recent heatwave
last_event_dates <- region_hw %>% 
  filter(mhw_event_no == last_event) %>% 
  pull(time)


# Buffer the dates, start 7 days before
event_start <- (min(last_event_dates) - 7)
event_stop  <- max(last_event_dates)
date_seq <- seq.Date(from = event_start,
                     to   = event_stop,
                     by   = 1)



# Load the heatwave dates
data_window <- data.frame(
  time = c(min(date_seq) , max(date_seq) ),
  lon  = crop_x,
  lat  = crop_y)


# Pull data off box
hw_stack <- oisst_window_load(data_window = data_window, 
                              anomalies = T, mac_os = "mojave")


#drop any empty years that bug in
hw_stack <- hw_stack[map(hw_stack, class) != "character"]


##### Format the layers and loop through the maps  ####

# Grab only current year, format dates
this_yr   <- stack(hw_stack)
day_count <- length(names(this_yr))
day_labs  <- str_replace_all(names(this_yr), "[.]","-")
day_labs  <- str_replace_all(day_labs, "X", "")
day_count <- c(1:day_count) %>% setNames(day_labs)

# Progress through daily timeline to indicate heatwave status and severity
hw_timeline <- region_hw %>% 
  filter(time %in% as.Date(day_labs))
####  Plot Settings:

# Set palette limits to center it on 0 with scale_fill_distiller
limit <- c(max(values(this_yr), na.rm = T) * -1, 
           max(values(this_yr), na.rm = T) )


# Plot Heatwave 1 day at a time as a GIF
day_plots <- imap(day_count, function(date_index, date_label) {
  
  # grab dates
  heatwaves_st  <- st_as_stars(this_yr[[date_index]])
  
  #### 1. Map the Anomalies in Space
  day_plot <- ggplot() +
    geom_stars(data = heatwaves_st) +
    geom_sf(data = new_england, fill = "gray30", color = "white", size = .25) +
    geom_sf(data = canada, fill = "gray30", color = "white", size = .25) +
    geom_sf(data = greenland, fill = "gray30", color = "white", size = .25) +
    geom_sf(data = region_extent, 
            color = gmri_cols("gmri blue"), 
            linetype = 2, size = 1,
            fill = "transparent") +
    scale_fill_distiller(palette = "RdYlBu", 
                         na.value = "transparent", 
                         limit = limit) +
    map_theme(legend.position = "bottom") +
    coord_sf(xlim = crop_x, 
             ylim = crop_y, 
             expand = T) +
    guides("fill" = guide_colorbar(
      title = "Sea Surface Temperature Anomaly \u00b0C",
      title.position = "top",
      title.hjust = 0.5,
      barwidth = unit(3, "in"),
      frame.colour = "black",
      ticks.colour = "black")) 
  
  # Set colors by name
  color_vals <- c(
    "Sea Surface Temperature" = "royalblue",
    "Heatwave Event"          = "darkred",
    "MHW Threshold"           = "coral3",
    "Daily Climatology"        = "gray30")
  
  
  
  #### 2.  Plot the day and the overall anomaly to track dates
  date_timeline <- ggplot(data = hw_timeline, aes(x = time)) +
    geom_line(aes(y = sst, color = "Sea Surface Temperature")) +
    geom_line(aes(y = hwe, color = "Heatwave Event")) +
    geom_line(aes(y = cse, color = "Cold Spell Event")) +
    geom_line(aes(y = mhw_thresh, color = "MHW Threshold"), lty = 3, size = .5) +
    geom_line(aes(y = mcs_thresh, color = "MCS Threshold"), lty = 3, size = .5) +
    geom_line(aes(y = seas, color = "Daily Climatology"), lty = 2, size = .5) +
    scale_color_manual(values = color_vals) +
    
    # Animated Point /  line
    geom_point(
      data = filter(hw_timeline, time == as.Date(date_label)),
      aes(time, sst, shape = factor(mhw_event)), 
      color = gmri_cols("gmri blue"), 
      size = 3, show.legend = FALSE) + 
    geom_vline(data = filter(hw_timeline, time == as.Date(date_label)),
               aes(xintercept = time), 
               color = "gray50", 
               size = 0.5,
               linetype = 3,
               alpha = 0.8) + 
    labs(x = "", 
         y = "",
         color = "",
         subtitle = "Regional Temperature \u00b0C",
         shape = "Heatwave Event") +
    theme(legend.position = "bottom")
  
  
  ####  3. Assemble plot(s)
  p_layout <- c(
    area(t = 1, l = 1, b = 2, r = 8),
    area(t = 3, l = 1, b = 8, r = 8))
  
  # plot_agg <- (date_timeline / day_plot) + plot_layout(heights = c(1, 3))
  plot_agg <- date_timeline + day_plot + plot_layout(design = p_layout)
  
  
  return(plot_agg )
  
  
})


walk(day_plots, print)

Ranking Warming Rates

If we look at the rates of change from 1982-2021 for each grid cell, rather than the observed temperature, it is possible to rank how hot each location on earth is warming relative to all the others.

Once we have the rankings, we can then take the average ranking within the Gulf of Maine we can obtain the average warming rank for the area compared to the rest of the globe.

# 1. Warming Rates and Rankings
rates_path <- paste0(oisst_path, "warming_rates/annual_warming_rates")
rates_stack_all <- stack(str_c(rates_path, "1982to2021.nc"), 
                         varname = "annual_warming_rate")
ranks_stack_all <- stack(str_c(rates_path, "1982to2021.nc"), 
                         varname = "rate_percentile")
# Get the rank information that go with the original extent used by the timelines
ranks_masked <- mask_nc(ranks_stack_all, region_extent)
rates_masked <- mask_nc(rates_stack_all, region_extent)
region_ranks <- get_masked_vals(ranks_masked, rates_masked)

# Prep it for text input.
avg_rank <- region_ranks$`Mean Rank` *100
avg_rate <- region_ranks$`Mean Rate`
low_rank <- region_ranks$`Min Rank` *100
low_rate <- region_ranks$`Min Rate`
top_rank <- region_ranks$`Max Rank` *100
top_rank <- ifelse(top_rank == 100, "greater than or equal to 99.5", top_rank)
top_rate <- region_ranks$`Max Rate`

Based on data from 1982-2021, the warming rates of Gulf Of Maine have been some of the highest in the world. The area as a whole has been increasing at a rate of 0.044\(^{\circ}C/year\) which is faster than 95.9% of the world’s oceans.

Over that same period locations within the Gulf Of Maine have been warming at rates as low as 0.017\(^{\circ}C/year\) and as rapidly as 0.094\(^{\circ}C/year\), corresponding to ranks as low as 60.3% and as high as greater than or equal to 99.5%.

Mapped below are the corresponding warming rates and their global rankings.

# For the full map  we mask again, but zoom out a little

# Mask again using the expanded mask so we can zoom out
ranks_masked <- mask_nc(ranks_stack_all, region_extent_expanded)
rates_masked <- mask_nc(rates_stack_all, region_extent_expanded)


# Make stars object
rank_stars <- st_as_stars(ranks_masked)
rates_stars <- st_as_stars(rates_masked)
####  Rates Map  ####

# Make Contours if desired
# rates_contour <-  st_contour(x = rates_stars, na.rm = T, breaks = seq(0.85, 1, by = 0.02))

# # rates map Original
# rates_map <- ggplot() +
#   geom_stars(data = rates_stars) +
#   geom_sf(data = new_england, fill = "gray90", size = .25) +
#   geom_sf(data = canada, fill = "gray90", size = .25) +
#   geom_sf(data = greenland, fill = "gray90", size = .25) +
#   geom_sf(data = region_extent, 
#           color = "black", 
#           fill = "transparent", linetype = 2, size = 0.5) +
#   scale_fill_viridis_c(option = "plasma", na.value = "transparent") +
#   map_theme +
#   coord_sf(xlim = crop_x, 
#            ylim = crop_y, expand = F) +
#   guides("fill" = guide_colorbar(
#     title = "Annual Temperature Change \u00b0C/year",
#     title.position = "top", 
#     title.hjust = 0.5,
#     barwidth = unit(2.5, "in"),
#     frame.colour = "black",
#     ticks.colour = "black"))



# Look at the spread
# hist(rates_stars$Annual.Sea.Surface.Temperature.Warming.Rate)

# Reclassify to discrete bins
rates_reclass <- rates_stars %>% 
  mutate(
    warm_rate = Annual.Sea.Surface.Temperature.Warming.Rate,
    rates_reclass = case_when(
      warm_rate > .09  ~ "> 0.1", 
      warm_rate > .08 ~ "0.08 - 0.1", 
      warm_rate > .06 ~ "0.06 - 0.08", 
      warm_rate > .04 ~ "0.05 - 0.06", 
      warm_rate > .04 ~ "0.04 - 0.05", 
      warm_rate > .02 ~ "0.02 - 0.04", 
      warm_rate < 0.2 ~ "< 0.02",
      TRUE ~ "NA"),
    rates_reclass = factor(rates_reclass, levels = c(
      "> 0.1", "0.08 - 0.1", 
      "0.06 - 0.08", "0.05 - 0.06", "0.04 - 0.05", 
      "0.02 - 0.04", "< 0.02",
      "Outside Area Measured"))
  ) %>% select(rates_reclass)


# Testing discrete scales
rates_map <- ggplot() +
  geom_stars(data = rates_reclass) +
  geom_sf(data = new_england, fill = "gray90", size = .25) +
  geom_sf(data = canada, fill = "gray90", size = .25) +
  geom_sf(data = greenland, fill = "gray90", size = .25) +
  geom_sf(data = region_extent, 
          color = "black", 
          fill = "transparent", linetype = 2, size = 0.5) +
  scale_fill_brewer(palette = "YlOrRd", 
                    na.value = "transparent", 
                    na.translate = F,
                    direction = -1) +
  coord_sf(xlim = crop_x, 
           ylim = crop_y, expand = F) +
  map_theme(legend.position = "bottom") +
  guides("fill" = guide_legend(
    title = "Annual Temperature Change \u00b0C/year",
    title.position = "top", 
    title.hjust = 0.5,
    label.position = "bottom", 
    keywidth = unit(1.5, "cm"),
    reverse = T,
    nrow = 1))


rates_map
grid::grid.raster(lab_logo, x = 0.05, y = 0.03, just = c('left', 'bottom'), width = unit(0.8, 'inches'))

####  Warming Ranks Map  ####

# Make Contours if desired
ranks_contour <-  st_contour(x = rank_stars, na.rm = T, breaks = seq(0.85, 1, by = 0.02))


# # Original Map
# # ranks map
# ranks_map <- ggplot() +
#   geom_stars(data = rank_stars) +
#   geom_sf(data = ranks_contour, fill = "transparent", color = "gray30", size = 0.1) +
#   geom_sf(data = new_england, fill = "gray90", size = .25) +
#   geom_sf(data = canada, fill = "gray90", size = .25) +
#   geom_sf(data = greenland, fill = "gray90", size = .25) +
#   geom_sf(data = region_extent, 
#           color = "black", 
#           fill = "transparent", linetype = 2, size = 0.5) +
#   scale_fill_viridis_c(option = "plasma",
#                        na.value = "transparent",
#                        limit = c(0.85, 1),
#                        oob = scales::oob_squish) +
#  map_theme +
#   coord_sf(xlim = crop_x, 
#            ylim = crop_y, expand = F) +
#   guides("fill" = guide_colorbar(
#     title = "Global Percentile of Warming Rates",
#     title.position = "top", 
#     title.hjust = 0.5,
#     barwidth = unit(2.5, "in"),
#     frame.colour = "black",
#     ticks.colour = "black")) +
#   labs(caption = "Ranking color scale truncated to display ranges of 0.85-1 
#                   Lower values will display as 0.85 or 85%
#                   Contour values every 2%")


# Look at the spread
# hist(rank_stars$Annual.Sea.Surface.Temperature.Warming.Rate.Rank)

# reclassify to discrete bins
ranks_reclass <- rank_stars %>% 
  mutate(
    warm_rank = Annual.Sea.Surface.Temperature.Warming.Rate.Rank,
    ranks_class = case_when(
      warm_rank > .999  ~ "> 99.9%", 
      warm_rank > .99  ~ "99 - 99.9%", 
      warm_rank > .98 ~ "98 - 99%", 
      warm_rank > .95 ~ "95 - 98%", 
      warm_rank > .90 ~ "90 - 95%", 
      warm_rank > .80 ~ "80 - 90%", 
      warm_rank < .80 ~ "< 80%",
      TRUE ~ "NA"),
    ranks_class = factor(ranks_class, levels = c(
      "> 99.9%", "99 - 99.9%", "98 - 99%", 
      "95 - 98%", "90 - 95%", 
      "80 - 90%", "< 80%"))
  ) %>% select(ranks_class)


ranks_map <- ggplot() +
  geom_stars(data = ranks_reclass) +
  geom_sf(data = new_england, fill = "gray90", size = .25) +
  geom_sf(data = canada, fill = "gray90", size = .25) +
  geom_sf(data = greenland, fill = "gray90", size = .25) +
  geom_sf(data = region_extent, 
          color = "black", 
          fill = "transparent", linetype = 2, size = 0.5) +
 scale_fill_brewer(palette = "YlOrRd", 
                    na.value = "transparent", 
                    na.translate = F,
                    direction = -1) +
  map_theme(legend.position = "bottom") +
  coord_sf(xlim = crop_x, 
           ylim = crop_y, expand = F) +
  guides("fill" = guide_legend(
    title = "Global Percentile of Warming Rates",
    title.position = "top", 
    title.hjust = 0.5,
    label.position = "bottom", 
    keywidth = unit(1.5, "cm"),
    reverse = T,
    nrow = 1)) +
  labs(caption = "Pixels ranked globally based on annual warming rate.")

# plot both
ranks_map
grid::grid.raster(lab_logo, x = 0.05, y = 0.03, just = c('left', 'bottom'), width = unit(0.8, 'inches'))

Shifting Baselines

In 2021 NOAA is transitioning standard climatologies from the 30-year period of 1982-2011 to a new period spanning 1992-2020. Changes in climate regimes often does not result in a uniform upward or downward change that is consistent throughout the year.

The plot below shows just how both the average temperature, as well as the annual highs and lows have shifted. When looking specifically at Gulf Of Maine here is how the expected temperature for each day of the year has shifted.

From this we can see that the Fall temperatures have risen more than those of the spring. There is also a large change in where the threshold for Marine Heatwave events sits, a consequence of exceptionally warm Fall temperatures becoming more common.

# Run heatwave detection using new climate period
heatwaves_91 <- pull_heatwave_events(region_timeseries, 
                                     threshold = 90, 
                                     clim_ref_period = c("1991-01-01", "2020-12-31")) 



# Subtract old from the new
heatwaves_91 <- heatwaves_91 %>% 
  mutate(clim_shift = seas - region_hw$seas,
         upper_shift = mhw_thresh - region_hw$mhw_thresh,
         lower_shift = mcs_thresh - region_hw$mcs_thresh)


# Make arrows where we want to point at things:
arrows <- 
  tibble(
    x1 = as.Date(c("2000-07-15")),
    x2 = as.Date(c("2000-08-28")),
    y1 = c(1.25), 
    y2 = c(1.15)
  )

# arrows



# Plot the differences
heatwaves_91 %>% 
  filter(time >= last_year) %>% 
  mutate(year = year(time),
         yday = yday(time),
         flat_date = as.Date(yday-1, origin = base_date)) %>% 
  distinct(flat_date, .keep_all = T) %>%        
  ggplot(aes(x = flat_date)) +
  geom_line(aes(y = clim_shift, color = "Mean Temperature Shift")) + 
  geom_line(aes(y = upper_shift, color = "MHW Threshold Change")) + 
  geom_line(aes(y = lower_shift, color = "MCS Threshold Change")) + 
  geom_curve(
    data = arrows, aes(x = x1, y = y1, xend = x2, yend = y2),
    arrow = arrow(length = unit(0.08, "inch")), size = 0.5,
    color = "gray20", curvature = -0.3) +
  annotate("text", x = as.Date("2000-06-01"), y = 1.225, label = "Fall Extremes\nMore Frequent") +
  labs(x = "", 
       y = "Shift in Expected Temperature \u00b0C",
       color = "") + 
  theme(legend.position = "bottom") +
  scale_color_gmri() +
  scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0))

 

A work by Adam A. Kemberling

Akemberling@gmri.org